%----------------------------
% clear plots
 close all;
%---------------------------

// Christoph Gortz and John Tsoukalas
// News and Financial Intermediation in Aggregate Fluctuations, Review of Economics and Statistics

// Two sector model with price and wage rigidities with permanent and transitory 
// technology shocks and separate specification of Multi Factor Productivity in I anc C sectors
// Capital is immobile across sectors and there are two investment and capital utilization decisions
// intratemporal investment adjustment cocts are incorporated a-la Huffman and Whynne
// This version does not include an Investment-specific shock v_l and a TFP transitory shock to cons sector z_l!
// This model includes a banking sector a la Gertler and Karadi (2010) 

var   labi labc lab kc ki kpi kpc mcc mci zcapi zcapc rki rkc pi lam phifi phifc c invei invec inve y pinfc pinfi w r 
z v b g qs ms sw ewma z_l v_l dy dc dinve    dpi   dw labobs robs pinfobsc pinfobsi spreadobsc spreadobsi rrate //spread dcredit
epinfmac epinfmai spinfc spinfi spi spkc spki
zcapif zcapcf rkcf rkif kpif kpcf kif kcf lamf phifif phifcf cf invef inveif invecf yf labf labif labcf pif wf 
sdfbf nucf nuif rbcf rbif z1cf z1if etaif etacf z2cf z2if lrcf lrif ncf nif qcf qif scf sif
sdfb rbc rbi z1c z1i z2c z2i nuc nui etac etai lrc lri qc qi sc si nc ni spreadc spreadi snc sni psic psii
epkc1aux epkc2aux epkc3aux epkc4aux epkc8aux epki1aux epki2aux epki3aux epki4aux epki8aux
ez1aux ez4aux ez8aux ev1aux ev4aux ev8aux 
cthetab dequity tfp etfp4aux etfp8aux rbcq rbiq
;


varexo ez ev eb em epinfc epinfi ew enc eni eg epkc epki eqs evl ezl espi 
epkc1 epkc2 epkc3 epkc4 epkc8 epki1 epki2 epki3 epki4 epki8 ez1 ez4 ez8 ev1 ev4 ev8 ethetab etfp etfp4 etfp8;
 
//esp1, esp2 are measurement errors

parameters constepinfc constepinfi constebeta cmaw cmapi cmapc  alfac alfai
czcapi czcapc czcap csadjcost cdeltai cdeltac csigma chabb  
cindw cprobw cindpc cprobpc cindpi cprobpi csigl cfc
crpi crdy cry crr 
crhoz crhov crhob crhog crhoqs crhoms crhopinf crhow crhozl crhovl crhospkc crhospki
cminrho
cg cga 
cgv clandapaux clandawaux
Lss nvalueic 
crhopinfc crhopinfi crhospi
cthetabss crhothetab clamb crhosnc crhosni sigmasn ctau cvarthetac cvarthetai psicss psiiss lridata lrcdata //equityss//cvaromega
grconstepinfc crhotfp csadjcostc csadjcosti

// Steady state parameters

X1aux crki crkc cpi cw KcLc KiLi Lc Li css iss Kc Ki K Kibar Kcbar yss iiss icss
// Parameters needed to calculate ss parameters
cbeta ga 
gv clandap clandaw
// Steady state parameters resulting from financial intermediaries
z1css z1iss z2css z2iss rbcss rbiss rss lrcss lriss nucss nuiss qcss qiss sncss sniss scss siss;


// fixed parameters

cdeltai =.025;   // depreciation rate I sector capital
cdeltac =.025;   // depreciation rate C sector capital
alfac=.3;
alfai=.3;

// SS restriction is cbeta = pinfobs/R, for US data (1990Q2:2011Q1) this gives cbeta = 0.9974
// from this you can get implied value for constebeta = (cbeta^-1 -1)*100
constebeta = 0.2607;   // equivalent to cbeta = 0.9974, constistent with 1990Q2-2011Q1 US data
constepinfc = 0.6722;  // constistent with 1990Q2-2011Q1 US data
grconstepinfc = constepinfc/100 + 1;
constepinfi = 0.0245;  // constistent with 1990Q2-2011Q1 US data
cg = 0.177;            // ss nominal government spending to nominal output ratio, constistent with 1982Q4-2011Q1 US data
nvalueic = 0.399;      // average value of nominal investment to consumption ratio, constistent with 1982Q4-2011Q1 US data
Lss = 1;               // ss of total hours (by choosing appropriate value of the weight of leisure in utility)
cmapc =    0.0;          // MA price mark-up
cmapi =    0.0;          // MA price mark-up
cmaw  =   0.0;           // MA wage mark-up
cminrho = 0.0;      // intrtemporal inv. adjustment cost, zero equals crho = -1
csigma=1.0;
cfc=1.5;




// estimated parameters initialisation

cga = 0.141; 
cgv = 0.434;

csadjcostc= 4;	// investment adj. cost parameter
csadjcosti= 4;	// investment adj. cost parameter

chabb=    0.5;      // consumption habit 
cprobw=   0.66;     // wage calvo contract
cindw=    0.5;      // wage indexation
csigl=    2;      // inverse labour supply ela
cprobpc=   0.66;    // consumption sector price calvo probability 
cindpc=    0.5;     // consumption sector price indexation
cprobpi=   0.66;    // investment sector price calvo probability 
cindpi=    0.5;     // investmnt sector price indexation
czcapi=    5;    // variable I-capital utilization 
czcapc=    5;    // variable C-capital utilization 
czcap =    5.75;
crpi=     1.7;      // taylor rule inflation 
crr=      0.6;        // taylor rule inertia  
cry=      0.0;        // taylor rule output gap 
crdy=     0.25;        // taylor rule output gap growth
clandawaux = 1.15-1; // ss wage mark up
clandapaux = 1.15-1; // ss price mark up

crhoz=    0.4;      // consumption sector productivity AR(1) --- growth shock
crhov=    0.4;        // investment sector productivity AR(1) --- growth shock
crhozl=    0.0;     // consumption sector productivity AR(1) --- stationary shock
crhovl=    0.0;     // investment sector productivity AR(1) --- stationary shock
crhob=    0.6;      // preference shock AR(1)
crhog=    0.6;      // gov spending shock AR(1)
crhoqs=   0.0;      // installation shock AR(1)
crhoms=   0.0;      // monetary policy shock AR(1)
crhopinfi= 0.6;     // price mark-up shock AR(1) in I sector
crhopinfc= 0.6;     // price mark-up shock AR(1) in C sector
crhow=    0.6;      // wage mark-up shock AR(1)
crhospkc = 0.6;     // quality of physical capital shock in C sector
crhospki = 0.6;     // quality of physical capital shock in I sector
crhospi = 0.0;      // persistence of inflation objective shock
crhosnc = 0.0;      // Persistence of shock to intermediarys equity capital in C-sector
crhosni = 0.0;      // Persistence of shock to intermediarys equity capital in I-sector
crhothetab = 0;     // Persistence of equity shock
crhotfp = 0.0;      // Persistence common TFP


// Fixed parameters specific to financial intermediaries
cthetabss = 0.96;      // Fraction of bankers that survive at following period (in SS)
lridata = 5.47;         // Leverage ratio implied by the data for I sector
lrcdata = lridata;      // Leverage ratio implied by the data for C sector
cvarthetac = 0;         // Response of gov's credit to external spread in C sector
cvarthetai = cvarthetac;// Response of gov's credit to external spread in I sector
psicss = 0;             // Steady state fraction assets intermed by CB in C sector
psiiss = psicss;        // Steady state fraction assets intermed by CB in I sector
ctau = 0.0;             // Government's cost of providing profit
//equityss = 0.5;         // C-sector equity over total equity in steady state
// Central bank as an intermediary is shut down for psicss = psiss = ctau = cvarthetac = cvarthetai = 0


// Steady state parameters
// Initial values according to steady state with rho = -1
// parameters needed for this section:
cbeta = 1/(1+constebeta/100);
ga = cga/100 ;
gv = cgv/100 ;
clandap = 1 + clandapaux;
clandaw = 1 + clandawaux;
//constepinfi  = constepinfc + (cga+ (alfac-1)/(1-alfai)*cgv);





model (linear); 

// Auxiliary parameter
//# crho = -(cminrho + 1);  // rho = -rho
# crho = 1/(cminrho - 1);   // rho = -rho
//# cga = constepinfi - constepinfc -  (alfac-1)/(1-alfai)*cgv;
//# ga = cga/100;



// flexible economy

          
          invef = clandap*(alfai*kif + (1-alfai)*labif + v_l);
          invef = (iiss^-crho + icss^-crho)^-1*iiss^-crho*inveif + (iiss^-crho + icss^-crho)^-1*icss^-crho*invecf;          
          cf + crki*Kibar/css*exp(-(1/(1-alfai))*gv)*zcapif + crkc*Kcbar/css*exp(-(1/(1-alfai))*gv)*zcapcf = clandap*(alfac*kcf + (1-alfac)*labcf + z_l);      
          z_l - alfac*labcf + alfac*kcf = wf;
          z_l + (1 - alfac)*labcf + (alfac-1)*kcf = rkcf;
          v_l - alfai*labif + alfai*kif + pif = wf;
          v_l + (1 - alfai)*labif + (alfai-1)*kif + pif = rkif;
          
      lamf = (chabb*cbeta*exp(ga+(alfac/(1-alfai))*gv)/((exp(ga+(alfac/(1-alfai))*gv)-chabb*cbeta)*(exp(ga+(alfac/(1-alfai))*gv)-chabb)))*cf(+1)       
        - ((exp(2*(ga+(alfac/(1-alfai))*gv)) + chabb^2*cbeta)/((exp(ga+(alfac/(1-alfai))*gv)-chabb*cbeta)*(exp(ga+(alfac/(1-alfai))*gv)-chabb)))*cf       
        + (chabb*exp(ga+(alfac/(1-alfai))*gv)/((exp(ga+(alfac/(1-alfai))*gv)-chabb*cbeta)*(exp(ga+(alfac/(1-alfai))*gv)-chabb)))*cf(-1)
        + ( exp(ga+(alfac/(1-alfai))*gv)*chabb*cbeta*crhoz - chabb*exp(ga+(alfac/(1-alfai))*gv) )/((exp(ga+(alfac/(1-alfai))*gv)-chabb*cbeta)*(exp(ga+(alfac/(1-alfai))*gv)-chabb))*z
        + b
        +  (alfac/(1-alfai))*( exp(ga+(alfac/(1-alfai))*gv)*chabb*cbeta*crhov - chabb*exp(ga+(alfac/(1-alfai))*gv) )/((exp(ga+(alfac/(1-alfai))*gv)-chabb*cbeta)*(exp(ga+(alfac/(1-alfai))*gv)-chabb))*v;

        // normalization: ( exp(ga+(alfac/(1-alfai))*gv) - chabb*cbeta*crhob )/(exp(ga+(alfac/(1-alfai))*gv)-chabb*cbeta)*b
        
        //phifif = (1-cdeltai)*cbeta*exp(-(1/(1-alfai))*gv)*( phifif(+1) - z(+1) - 1/(1-alfai)*v(+1) )
        //+ (1 - (1-cdeltai)*cbeta*exp(-(1/(1-alfai))*gv))*( lamf(+1) - z(+1) - 1/(1-alfai)*v(+1) + rkif(+1) );

        //phifcf = (1-cdeltac)*cbeta*exp(-(1/(1-alfai))*gv)*( phifcf(+1) - z(+1) - 1/(1-alfai)*v(+1) )
        //+ (1 - (1-cdeltac)*cbeta*exp(-(1/(1-alfai))*gv))*( lamf(+1) - z(+1) - 1/(1-alfai)*v(+1) + rkcf(+1) );
        
        rkif = czcapi*zcapif;

        rkcf = czcapc*zcapcf;
        
       
        lamf = -pif + phifif + qs - exp(2*((1/(1-alfai))*gv))*csadjcostc*(inveif - inveif(-1) + z + 1/(1-alfai)*v) 
        + cbeta*exp(2*((1/(1-alfai))*gv))*csadjcostc*( inveif(+1) - inveif + z(+1) +  1/(1-alfai)*v(+1) ) -(1+crho)*((iiss^-crho + icss^-crho)^-1*(icss^-crho*invecf+iiss^-crho*inveif)-inveif);

        lamf = -pif + phifcf + qs - exp(2*((1/(1-alfai))*gv))*csadjcostc*(invecf - invecf(-1) + z + 1/(1-alfai)*v) 
        + cbeta*exp(2*((1/(1-alfai))*gv))*csadjcostc*( invecf(+1) - invecf + z(+1) +  1/(1-alfai)*v(+1) ) -(1+crho)*((iiss^-crho + icss^-crho)^-1*(icss^-crho*invecf+iiss^-crho*inveif)-invecf);
         
        
	    kif =  kpif(-1)+ zcapif + spki  - 1/(1-alfai)*v;

        kcf =  kpcf(-1)+ zcapcf + spkc  - 1/(1-alfai)*v;
        
        kpif =  (1-cdeltai)*exp(-(1/(1-alfai))*gv)*( kpif(-1) + spki  - (1/(1-alfai))*v ) + ( 1 - (1-cdeltai)*exp(-(1/(1-alfai))*gv))*(qs + inveif);

        kpcf =  (1-cdeltac)*exp(-(1/(1-alfai))*gv)*( kpcf(-1) + spkc  - (1/(1-alfai))*v ) + ( 1 - (1-cdeltac)*exp(-(1/(1-alfai))*gv))*(qs + invecf);

        wf =  b + csigl*labf - lamf;

        yf = css/yss*cf+cpi*iss/yss*(invef + pif) + g ;
        //(1-cg)*yf = (css/yss)*cf+(cpi*iss/yss)*(invef+ pif) + (1-cg)*g + ctau/yss* ( psicss*qcss*scss*(psic+qcf+scf) + psiiss*qiss*siss*(psii+qif+sif) ); 

        labf = Lc/Lss*labcf + Li/Lss*labif;
        

        //Equations needed to implement financial intermediaries

        // Stochastic discount factor of bank
        sdfbf = lamf-lamf(-1);

        // Definition of nu for C sector
		//nucf = (1-cthetab*cbeta*z1css)*(sdfbf(+1)-z(+1)-alfac/(1-alfai)*v(+1)) + (1-cthetab*cbeta*z1css)/(rbcss-rss)*(rbcss*rbcf(+1) - rss*r)  +  cthetab*cbeta*z1css*(z1cf(+1)+nucf(+1));
        //NEW - with stochastic cthetab:
        nucf = -cthetabss*(1-cthetabss*cbeta*z1css)/(1-cthetabss) * cthetab + (1-cthetabss*cbeta*z1css)*(sdfbf(+1)-z(+1)-alfac/(1-alfai)*v(+1)) + (1-cthetabss*cbeta*z1css)/(rbcss-rss)*(rbcss*rbcf(+1) - rss*r)  +  cthetabss*cbeta*z1css*(z1cf(+1)+nucf(+1)+cthetab);
        
        // Definition of nu for I sector
		//nuif = (1-cthetab*cbeta*z1iss)*(sdfbf(+1)-z(+1)-alfac/(1-alfai)*v(+1)) + (1-cthetab*cbeta*z1iss)/(rbiss-rss)*(rbiss*rbif(+1) - rss*r)  +  cthetab*cbeta*z1iss*(z1if(+1)+nuif(+1));
        //NEW - with stochastic cthetab: 
        nuif = -cthetabss*(1-cthetabss*cbeta*z1iss)/(1-cthetabss) * cthetab + (1-cthetabss*cbeta*z1iss)*(sdfbf(+1)-z(+1)-alfac/(1-alfai)*v(+1)) + (1-cthetabss*cbeta*z1iss)/(rbiss-rss)*(rbiss*rbif(+1) - rss*r)  +  cthetabss*cbeta*z1iss*(z1if(+1)+nuif(+1)+cthetab);
        
        // Definition of eta for C sector
		//etacf = (1-cthetab*cbeta*z2css)*(sdfbf(+1)-z(+1)-alfac/(1-alfai)*v(+1) + r)  +   cthetab*cbeta*z2css*(z2cf(+1)+etacf(+1));
        //NEW - with stochastic cthetab:
        etacf = -cthetabss*(1-cthetabss*cbeta*z2css)/(1-cthetabss) * cthetab + (1-cthetabss*cbeta*z2css)*(sdfbf(+1)-z(+1)-alfac/(1-alfai)*v(+1) + r)  +   cthetabss*cbeta*z2css*(z2cf(+1)+etacf(+1)+cthetab);

        // Definition of eta for I sector
		//etaif = (1-cthetab*cbeta*z2iss)*(sdfbf(+1)-z(+1)-alfac/(1-alfai)*v(+1) + r)  +   cthetab*cbeta*z2iss*(z2if(+1)+etaif(+1));
        //NEW - with stochastic cthetab: 
        etaif = -cthetabss*(1-cthetabss*cbeta*z2iss)/(1-cthetabss) * cthetab + (1-cthetabss*cbeta*z2iss)*(sdfbf(+1)-z(+1)-alfac/(1-alfai)*v(+1) + r)  +   cthetabss*cbeta*z2iss*(z2if(+1)+etaif(+1)+cthetab);


        ////////////////////////////////
        // Alternative definitions for eta and nu according to Ferman's code
        //nucf = ((1-cbeta*1*exp(-ga-alfac/(1-alfai)*gv)*cthetab*z1css) + (  cthetab*z1css ))  *  (sdfbf(+1)-z(+1)-alfac/(1-alfai)*v(+1))   +   (1-cbeta*1*exp(-ga-alfac/(1-alfai)*gv)*cthetab*z1css)/(rbcss-rss) * (rbcss*rbcf(+1)-rss*r)    +    cbeta*1*exp(-ga-alfac/(1-alfai)*gv)*cthetab*z1css* (z1cf(+1)+nucf(+1));
        //nuif = ((1-cbeta*1*exp(-ga-alfac/(1-alfai)*gv)*cthetab*z1iss) + (  cthetab*z1iss ))  *  (sdfbf(+1)-z(+1)-alfac/(1-alfai)*v(+1))   +   (1-cbeta*1*exp(-ga-alfac/(1-alfai)*gv)*cthetab*z1iss)/(rbiss-rss) * (rbiss*rbif(+1)-rss*r)    +    cbeta*1*exp(-ga-alfac/(1-alfai)*gv)*cthetab*z1iss* (z1if(+1)+nuif(+1));
        //etacf = cthetab*cbeta*1*exp(-ga-alfac/(1-alfai)*gv)*z2css  *  (sdfbf(+1)-z(+1)-alfac/(1-alfai)*v(+1)+z2cf(+1)+etacf(+1) );
        //etaif = cthetab*cbeta*1*exp(-ga-alfac/(1-alfai)*gv)*z2iss  *  (sdfbf(+1)-z(+1)-alfac/(1-alfai)*v(+1)+z2if(+1)+etaif(+1) );
        // Alternative definnitions for eta and nu according to Gertler and Karadi's paper
        //nucf = (1-cthetab*cbeta*z1css*exp(-ga-alfac/(1-alfai)*gv))*(sdfbf(+1)-z(+1)-alfac/(1-alfai)*v(+1)) + (1-cthetab*cbeta*z1css*exp(-ga-alfac/(1-alfai)*gv))/(rbcss-rss)*(rbcss*rbcf(+1) - rss*r)  +  cthetab*exp(-ga-alfac/(1-alfai)*gv)*cbeta*z1css*(sdfbf(+1)-z(+1)-alfac/(1-alfai)*v(+1)+z1cf(+1)+nucf(+1));
        //nuif = (1-cthetab*cbeta*z1iss*exp(-ga-alfac/(1-alfai)*gv))*(sdfbf(+1)-z(+1)-alfac/(1-alfai)*v(+1)) + (1-cthetab*cbeta*z1css*exp(-ga-alfac/(1-alfai)*gv))/(rbiss-rss)*(rbiss*rbif(+1) - rss*r)  +  cthetab*exp(-ga-alfac/(1-alfai)*gv)*cbeta*z1iss*(sdfbf(+1)-z(+1)-alfac/(1-alfai)*v(+1)+z1if(+1)+nuif(+1));
        //etacf = (cthetab*cbeta*z2css*exp(-ga-alfac/(1-alfai)*gv))*(sdfbf(+1)-z(+1)-alfac/(1-alfai)*v(+1) +z2cf(+1)+etacf(+1) ) ;
        //etaif = (cthetab*cbeta*z2iss*exp(-ga-alfac/(1-alfai)*gv))*(sdfbf(+1)-z(+1)-alfac/(1-alfai)*v(+1) +z2if(+1)+etaif(+1) ) ;
        ////////////////////////////////////


        // Definition of z1 for C sector
		z1cf = lrcf - lrcf(-1) + z2cf;

        // Definition of z1 for I sector
		z1if = lrif - lrif(-1) + z2if;

        // Definition of z2 for C sector
		z2cf = 1/((rbcss-rss)*lrcss+rss) * ( rbcss*lrcss*rbcf + rss*(1-lrcss)*r(-1) + (rbcss-rss)*lrcss*lrcf(-1) );

        // Definition of z2 for I sector
		z2if = 1/((rbiss-rss)*lriss+rss) * ( rbiss*lriss*rbif + rss*(1-lriss)*r(-1) + (rbiss-rss)*lriss*lrif(-1) );

        // Leverage ratio of bank in C sector
		lrcf = etacf + nucss/(clamb-nucss)*nucf;

        // Leverage ratio of bank in I sector
		lrif = etaif + nuiss/(clamb-nuiss)*nuif;

        // Wealth accumulation equation of bank in C sector
    	//ncf = sncss*cthetab*lrcss*exp(-ga-alfac/(1-alfai)*gv) * (rbcss*rbcf + (1/lrcss-1)*rss*r(-1) + (rbcss-rss)*lrcf(-1))    +      sncss*cthetab*exp(-ga-alfac/(1-alfai)*gv) * (  (rbcss-rss)*lrcss + rss  ) * (  -z-alfac/(1-alfai)*v+ncf(-1)  )     +     (1-sncss*cthetab*exp(-ga-alfac/(1-alfai)*gv) * (  (rbcss-rss)*lrcss + rss ) ) * (qcf+scf - psicss/(1-psicss)*psic )  +  (  cthetab*exp(-ga-alfac/(1-alfai)*gv)*((rbcss-rss)*lrcss+rss)  + ( 1-cthetab*((rbcss-rss)*lrcss+rss)  ) ) * snc;
        // NEW - with stochastic cthetab:
        ncf = sncss*cthetabss*exp(-ga-alfac/(1-alfai)*gv)*((rbcss-rss)*lrcss+rss) * (ncf(-1) + cthetab(-1) - z -alfac/(1-alfai)*v + snc)
            + sncss*cthetabss*exp(-ga-alfac/(1-alfai)*gv) * ( rbcss*lrcss*rbcf - rss*lrcss*r(-1) + (rbcss-rss)*lrcss*lrcf(-1) + rss*r(-1) )
            + (1-exp(-ga-alfac/(1-alfai)*gv)*((rbcss-rss)*lrcss+rss)*cthetabss*sncss) * (qcf + scf + snc);

        // Wealth accumulation equation of bank in I sector
    	//nif = sniss*cthetab*lriss*exp(-ga-alfac/(1-alfai)*gv) * (rbiss*rbif + (1/lriss-1)*rss*r(-1) + (rbiss-rss)*lrif(-1))    +      sniss*cthetab*exp(-ga-alfac/(1-alfai)*gv) * (  (rbiss-rss)*lriss + rss  ) * (  -z-alfac/(1-alfai)*v+nif(-1)  )     +     (1-sniss*cthetab*exp(-ga-alfac/(1-alfai)*gv) * (  (rbiss-rss)*lriss + rss ) ) * (qif+sif - psiiss/(1-psiiss)*psii )  +  (  cthetab*exp(-ga-alfac/(1-alfai)*gv)*((rbiss-rss)*lriss+rss)  + ( 1-cthetab*((rbiss-rss)*lriss+rss)  ) ) * sni;
        //NEW - with stochastic cthetab:
        nif = sniss*cthetabss*exp(-ga-alfac/(1-alfai)*gv)*((rbiss-rss)*lriss+rss) * (nif(-1) + cthetab(-1) - z -alfac/(1-alfai)*v + sni)
            + sniss*cthetabss*exp(-ga-alfac/(1-alfai)*gv) * ( rbiss*lriss*rbif - rss*lriss*r(-1) + (rbiss-rss)*lriss*lrif(-1) + rss*r(-1) )
            + (1-exp(-ga-alfac/(1-alfai)*gv)*((rbiss-rss)*lriss+rss)*cthetabss*sniss) * (qif + sif + sni);

        // Borrow in advance constraint for C sector
		kpcf = scf;

        // Borrow in advance constraint for I sector
		kpif = sif;

        // Stochastic return on assets for bank in C sector
    	rbcf = 1/(crkc+qcss*(1-cdeltac)) * ( crkc*(rkcf+zcapcf) + qcss*(1-cdeltac)*(qcf) ) -qcf(-1) +spkc + z - (1-alfac)/(1-alfai)*v;

        // Stochastic return on assets for bank in I sector
    	rbif = 1/(crki+qiss*(1-cdeltai)) * ( crki*(rkif+zcapif) + qiss*(1-cdeltai)*(qif) ) -qif(-1) +spki + z - (1-alfac)/(1-alfai)*v;

        // Price of capital in C sector
		qcf = exp(2*1/(1-alfai)*gv)*csadjcostc * ( invecf-invecf(-1)+1/(1-alfai)*v) - cbeta*exp(2*1/(1-alfai)*gv)*csadjcostc * (invecf(+1)-invecf+1/(1-alfai)*v(+1)) - qs + pif + (1+crho) * ( (iiss^(-crho)+icss^(-crho))^(-1) * (icss^(-crho)*invecf + iiss^(-crho)*inveif) - invecf );

        // Price of capital in I sector
		qif = exp(2*1/(1-alfai)*gv)*csadjcostc * ( inveif-inveif(-1)+1/(1-alfai)*v) - cbeta*exp(2*1/(1-alfai)*gv)*csadjcostc * (inveif(+1)-inveif+1/(1-alfai)*v(+1)) - qs + pif + (1+crho) * ( (iiss^(-crho)+icss^(-crho))^(-1) * (icss^(-crho)*invecf + iiss^(-crho)*inveif) - inveif );
        
        // Leverage equation in C sector
        qcf = lrcf + ncf - scf + psicss/(1-psicss)*psic;

        // Leverage equation in I sector
        qif = lrif + nif - sif + psiiss/(1-psiiss)*psii;

       

// sticky price - wage economy



          // consumption and investment sector output and cost minimization

          inve = clandap*(alfai*ki + (1-alfai)*labi + v_l + tfp);
          inve = (iiss^-crho + icss^-crho)^-1*iiss^-crho*invei + (iiss^-crho + icss^-crho)^-1*icss^-crho*invec;
          
          //inve = (iiss/iss)*invei + (icss/iss)*invec;
          c  + crki*Kibar/css*exp(-(1/(1-alfai))*gv)*zcapi + crkc*Kcbar/css*exp(-(1/(1-alfai))*gv)*zcapc = clandap*(alfac*kc + (1-alfac)*labc + z_l + tfp);

	      mcc =  alfac*rkc+(1-alfac)*w - z_l - tfp;
          mci =  alfai*rki+(1-alfai)*w - pi - v_l - tfp;
          rkc - w = labc - kc;
          rki - w = labi - ki;

          
        



         // consumption and investment sector optimal pricing --- Philips curves
          
           pinfc -spi =  (cbeta/(1+cindpc*cbeta))*(pinfc(+1)-spi) + (cindpc/(1+cindpc*cbeta))*(pinfc(-1)-spi) 
               +( (1-cprobpc)*(1-cbeta*cprobpc)/((1+cindpc*cbeta)*cprobpc) )* mcc  +  spinfc ;
               
           pinfi -spi =  (cbeta/(1+cindpi*cbeta))*(pinfi(+1)-spi) + (cindpi/(1+cindpi*cbeta))*(pinfi(-1)-spi) 
               +( (1-cprobpi)*(1-cbeta*cprobpi)/((1+cindpi*cbeta)*cprobpi) )* mci  + spinfi ; 

          
          // normalization of shocks:
          // spinf* = ( (1-cprobpi)*(1-cbeta*cprobpi)/((1+cindpi*cbeta)*cprobpi) )* spinf ; 
            
        // Households
        
        lam = (chabb*cbeta*exp(ga+(alfac/(1-alfai))*gv)/((exp(ga+(alfac/(1-alfai))*gv)-chabb*cbeta)*(exp(ga+(alfac/(1-alfai))*gv)-chabb)))*c(+1)
        - ((exp(2*(ga+(alfac/(1-alfai))*gv)) + chabb^2*cbeta)/((exp(ga+(alfac/(1-alfai))*gv)-chabb*cbeta)*(exp(ga+(alfac/(1-alfai))*gv)-chabb)))*c
        + (chabb*exp(ga+(alfac/(1-alfai))*gv)/((exp(ga+(alfac/(1-alfai))*gv)-chabb*cbeta)*(exp(ga+(alfac/(1-alfai))*gv)-chabb)))*c(-1)
        + ( exp(ga+(alfac/(1-alfai))*gv)*chabb*cbeta*crhoz - chabb*exp(ga+(alfac/(1-alfai))*gv) )/((exp(ga+(alfac/(1-alfai))*gv)-chabb*cbeta)*(exp(ga+(alfac/(1-alfai))*gv)-chabb))*z
	    + b
        +  (alfac/(1-alfai))*( exp(ga+(alfac/(1-alfai))*gv)*chabb*cbeta*crhov - chabb*exp(ga+(alfac/(1-alfai))*gv) )/((exp(ga+(alfac/(1-alfai))*gv)-chabb*cbeta)*(exp(ga+(alfac/(1-alfai))*gv)-chabb))*v;
        
	    //normalization of shocks:
        //b* =  ( exp(ga+(alfac/(1-alfai))*gv) - chabb*cbeta*crhob )/(exp(ga+(alfac/(1-alfai))*gv)-chabb*cbeta)*b

        lam = r + lam(+1) - z(+1) - v(+1)*alfac/(1-alfai) - pinfc(+1);
        
        //phifi = (1-cdeltai)*cbeta*exp(-(1/(1-alfai))*gv)*( phifi(+1) - 1/(1-alfai)*v(+1) )
        //+ (1 - (1-cdeltai)*cbeta*exp(-(1/(1-alfai))*gv))*( lam(+1) - 1/(1-alfai)*v(+1) + rki(+1) );

        //phifc = (1-cdeltac)*cbeta*exp(-(1/(1-alfai))*gv)*( phifc(+1) - 1/(1-alfai)*v(+1) )
        //+ (1 - (1-cdeltac)*cbeta*exp(-(1/(1-alfai))*gv))*( lam(+1) - 1/(1-alfai)*v(+1) + rkc(+1) );
        
        rki = czcapi*zcapi;
        rkc = czcapc*zcapc;

        //zcapi =  (1/(czcap/(1-czcap)))* rki ;
        //zcapc =  (1/(czcap/(1-czcap)))* rkc ;
        
        lam = -pi + phifi + qs - exp(2*((1/(1-alfai))*gv))*csadjcostc*(invei - invei(-1) + 1/(1-alfai)*v) 
        + cbeta*exp(2*((1/(1-alfai))*gv))*csadjcostc*( invei(+1) - invei +  1/(1-alfai)*v(+1) ) -(1+crho)*((iiss^-crho + icss^-crho)^-1* (icss^-crho*invec+iiss^-crho*invei)-invei);

        lam = -pi + phifc + qs - exp(2*((1/(1-alfai))*gv))*csadjcostc*(invec - invec(-1) + 1/(1-alfai)*v) 
        + cbeta*exp(2*((1/(1-alfai))*gv))*csadjcostc*( invec(+1) - invec +  1/(1-alfai)*v(+1) ) -(1+crho)*((iiss^-crho + icss^-crho)^-1* (icss^-crho*invec+iiss^-crho*invei)-invec);
        

       
	    ki =  kpi(-1)+ zcapi + spki - 1/(1-alfai)*v;

        kc =  kpc(-1)+ zcapc + spkc - 1/(1-alfai)*v;
 
        kpi =  (1-cdeltai)*exp(-(1/(1-alfai))*gv)*( kpi(-1) + spki - (1/(1-alfai))*v ) + ( 1 - (1-cdeltai)*exp(-(1/(1-alfai))*gv))*(qs + invei);

        kpc =  (1-cdeltac)*exp(-(1/(1-alfai))*gv)*( kpc(-1) + spkc - (1/(1-alfai))*v ) + ( 1 - (1-cdeltac)*exp(-(1/(1-alfai))*gv))*(qs + invec);
	       
	      w =  (1/(1+cbeta))*w(-1)
               +(cbeta/(1+cbeta))*w(+1)
               - ( (1-cprobw*cbeta)*(1-cprobw)/(cprobw*(1+cbeta)*(1+csigl*(1+1/(clandawaux)))))* (w - csigl*lab - b + lam)
               + (cindw/(1+cbeta))*(pinfc(-1)-spi)
               - ((1+cbeta*cindw)/(1+cbeta))*(pinfc-spi)
               + (cbeta/(1+cbeta))*(pinfc(+1)-spi)
               + (cindw/(1+cbeta))*( z(-1) + alfac/(1-alfai)*v(-1) )
               - (1+cbeta*cindw-crhoz*cbeta)/(1+cbeta)*z
               - ((1+cbeta*cindw-crhov*cbeta)/(1+cbeta))*(alfac/(1-alfai))*v
               + sw ;
               
      // normalization of shocks:
      //       sw* = ( (1-cprobw*cbeta)*(1-cprobw)/(cprobw*(1+cbeta)*(1+csigl*(1+1/(clandaw-1)))))* sw 
               
      lab = Lc/Lss*labc + Li/Lss*labi;
               
      y = (css/yss)*c+(cpi*iss/yss)*(inve+ pi) + g ;
      //(1-cg)*y = (css/yss)*c+(cpi*iss/yss)*(inve+ pi) + (1-cg)*g + ctau/yss* ( psicss*qcss*scss*(psic+qc+sc) + psiiss*qiss*siss*(psii+qi+si) );         
        
      // Taylor rule         
	  //    r =  crpi*(1-crr)*pinfc
      //         +cry*(1-crr)*(y-yf)     
      //         +crdy*(y-yf-y(-1)+yf(-1))
      //         +crr*r(-1)
      //        +ms  ;
    
    // ORIGINAL TAYLOR RULE

         // r =  crpi*(1-crr)*pinfc
         //      +cry*(1-crr)*(pinfc - pinfc(-1))
         //      +crdy*(1-crr)*(y-y(-1))
         //      +crr*r(-1)
         //     +ms  ;

        r =  crpi*(1-crr)*pinfc
               +cry*(1-crr)*0
               +crdy*(1-crr)*(y-y(-1))
               +crr*r(-1)
              +ms  ;

    //  Taylor rule similar to Smets Wouters 2005 to account for trend in interest rate observable
//        r =  spi + (1-crr)*( crpi*(pinfc(-1)-spi(-1)))
//               +(1-crr)*cry*(y-yf)     
//               +crdy*(y-yf-y(-1)+yf(-1))
//               +crr*(r(-1)-spi(-1))
//               +ms  ;
               
    pinfi - pinfc  = pi - pi(-1); 


   //dpi + (constepinfc - constepinfi) = pinfi  - pinfc;
// this implies the restriction
// constepinfi - constepinfc = (cga+ (alfac-1)/(1-alfai)*cgv)
// so initial parameter values need to respect this
// in order to solve the SS 

    //Equations needed to implement financial intermediaries
        // Stochastic discount factor of bank
        sdfb = lam-lam(-1);

        // Definition of nu for C sector
		//nuc = (1-cthetab*cbeta*z1css)*(sdfb(+1)-z(+1)-alfac/(1-alfai)*v(+1)) + (1-cthetab*cbeta*z1css)/(rbcss-rss)*(rbcss*rbc(+1) - rss*r)  +  cthetab*cbeta*z1css*(z1c(+1)+nuc(+1));
        //NEW - with stochastic cthetab:
        //nuc = -cthetabss*(1-cthetabss*cbeta*z1css)/(1-cthetabss) * cthetab + (1-cthetabss*cbeta*z1css)*(sdfb(+1)-z(+1)-alfac/(1-alfai)*v(+1)) + (1-cthetabss*cbeta*z1css)/(rbcss-rss)*(rbcss*rbc(+1) - rss*r)  +  cthetabss*cbeta*z1css*(z1c(+1)+nuc(+1)+cthetab);
        nuc = -cthetabss*(1-cthetabss*cbeta*z1css)/(1-cthetabss) * cthetab + (1-cthetabss*cbeta*z1css)*(sdfb(+1)-z(+1)-alfac/(1-alfai)*v(+1)) + (1-cthetabss*cbeta*z1css)/(rbcss*grconstepinfc-rss)*(rbcss*grconstepinfc*rbc(+1) + rbcss*grconstepinfc*pinfc(+1) - rss*r)  +  cthetabss*cbeta*z1css*(z1c(+1)+nuc(+1)+cthetab);
        
        // Definition of nu for I sector
		//nui = (1-cthetab*cbeta*z1iss)*(sdfb(+1)-z(+1)-alfac/(1-alfai)*v(+1)) + (1-cthetab*cbeta*z1iss)/(rbiss-rss)*(rbiss*rbi(+1) - rss*r)  +  cthetab*cbeta*z1iss*(z1i(+1)+nui(+1));
        //NEW - with stochastic cthetab: 
        //nui = -cthetabss*(1-cthetabss*cbeta*z1iss)/(1-cthetabss) * cthetab + (1-cthetabss*cbeta*z1iss)*(sdfb(+1)-z(+1)-alfac/(1-alfai)*v(+1)) + (1-cthetabss*cbeta*z1iss)/(rbiss-rss)*(rbiss*rbi(+1) - rss*r)  +  cthetabss*cbeta*z1iss*(z1i(+1)+nui(+1)+cthetab);
        nui = -cthetabss*(1-cthetabss*cbeta*z1iss)/(1-cthetabss) * cthetab + (1-cthetabss*cbeta*z1iss)*(sdfb(+1)-z(+1)-alfac/(1-alfai)*v(+1)) + (1-cthetabss*cbeta*z1iss)/(rbiss*grconstepinfc-rss)*(rbiss*grconstepinfc*rbi(+1) + rbiss*grconstepinfc*pinfc(+1) - rss*r)  +  cthetabss*cbeta*z1iss*(z1i(+1)+nui(+1)+cthetab);
        
        // Definition of eta for C sector
		//etac = (1-cthetab*cbeta*z2css)*(sdfb(+1)-z(+1)-alfac/(1-alfai)*v(+1) + r)  +   cthetab*cbeta*z2css*(z2c(+1)+etac(+1));
        //NEW - with stochastic cthetab:
        etac = -cthetabss*(1-cthetabss*cbeta*z2css)/(1-cthetabss) * cthetab + (1-cthetabss*cbeta*z2css)*(sdfb(+1)-z(+1)-alfac/(1-alfai)*v(+1) + r)  +   cthetabss*cbeta*z2css*(z2c(+1)+etac(+1)+cthetab);


        // Definition of eta for I sector
		//etai = (1-cthetab*cbeta*z2iss)*(sdfb(+1)-z(+1)-alfac/(1-alfai)*v(+1) + r)  +   cthetab*cbeta*z2iss*(z2i(+1)+etai(+1));
        //NEW - with stochastic cthetab: 
        etai = -cthetabss*(1-cthetabss*cbeta*z2iss)/(1-cthetabss) * cthetab + (1-cthetabss*cbeta*z2iss)*(sdfb(+1)-z(+1)-alfac/(1-alfai)*v(+1) + r)  +   cthetabss*cbeta*z2iss*(z2i(+1)+etai(+1)+cthetab);


        ////////////////////////////////
        // Alternative definitions for eta and nu according to Ferman's code
        //nuc = ((1-cbeta*1*exp(-ga-alfac/(1-alfai)*gv)*cthetab*z1css) + (  cthetab*z1css ))  *  (sdfb(+1)-z(+1)-alfac/(1-alfai)*v(+1))   +   (1-cbeta*1*exp(-ga-alfac/(1-alfai)*gv)*cthetab*z1css)/(rbcss-rss) * (rbcss*rbc(+1)-rss*r)    +    cbeta*1*exp(-ga-alfac/(1-alfai)*gv)*cthetab*z1css* (z1c(+1)+nuc(+1));
        //nui = ((1-cbeta*1*exp(-ga-alfac/(1-alfai)*gv)*cthetab*z1iss) + (  cthetab*z1iss ))  *  (sdfb(+1)-z(+1)-alfac/(1-alfai)*v(+1))   +   (1-cbeta*1*exp(-ga-alfac/(1-alfai)*gv)*cthetab*z1iss)/(rbiss-rss) * (rbiss*rbi(+1)-rss*r)    +    cbeta*1*exp(-ga-alfac/(1-alfai)*gv)*cthetab*z1iss* (z1i(+1)+nui(+1));
        //etac = cthetab*cbeta*1*exp(-ga-alfac/(1-alfai)*gv)*z2css  *  (sdfb(+1)-z(+1)-alfac/(1-alfai)*v(+1)+z2c(+1)+etac(+1) );
        //etai = cthetab*cbeta*1*exp(-ga-alfac/(1-alfai)*gv)*z2iss  *  (sdfb(+1)-z(+1)-alfac/(1-alfai)*v(+1)+z2i(+1)+etai(+1) );
        // Alternative definnitions for eta and nu according to Gertler and Karadi's paper
        //nuc = (1-cthetab*cbeta*z1css*exp(-ga-alfac/(1-alfai)*gv))*(sdfb(+1)-z(+1)-alfac/(1-alfai)*v(+1)) + (1-cthetab*cbeta*z1css*exp(-ga-alfac/(1-alfai)*gv))/(rbcss-rss)*(rbcss*rbc(+1) - rss*r)  +  cthetab*exp(-ga-alfac/(1-alfai)*gv)*cbeta*z1css*(sdfb(+1)-z(+1)-alfac/(1-alfai)*v(+1)+z1c(+1)+nuc(+1));
        //nui = (1-cthetab*cbeta*z1iss*exp(-ga-alfac/(1-alfai)*gv))*(sdfb(+1)-z(+1)-alfac/(1-alfai)*v(+1)) + (1-cthetab*cbeta*z1css*exp(-ga-alfac/(1-alfai)*gv))/(rbiss-rss)*(rbiss*rbi(+1) - rss*r)  +  cthetab*exp(-ga-alfac/(1-alfai)*gv)*cbeta*z1iss*(sdfb(+1)-z(+1)-alfac/(1-alfai)*v(+1)+z1i(+1)+nui(+1));
        //etac = (cthetab*cbeta*z2css*exp(-ga-alfac/(1-alfai)*gv))*(sdfb(+1)-z(+1)-alfac/(1-alfai)*v(+1) +z2c(+1)+etac(+1) ) ;
        //etai = (cthetab*cbeta*z2iss*exp(-ga-alfac/(1-alfai)*gv))*(sdfb(+1)-z(+1)-alfac/(1-alfai)*v(+1) +z2i(+1)+etai(+1) ) ;
        ////////////////////////////////////



        // Definition of z1 for C sector
		z1c = lrc - lrc(-1) + z2c;

        // Definition of z1 for I sector
		z1i = lri - lri(-1) + z2i;

        // Definition of z2 for C sector
		//z2c = 1/((rbcss-rss)*lrcss+rss) * ( rbcss*lrcss*rbc + rss*(1-lrcss)*r(-1) + (rbcss-rss)*lrcss*lrc(-1) );
        //z2c = grconstepinfc/((rbcss-rss)*lrcss+rss) * ( rbcss*lrcss*(rbc+pinfc) + rss/grconstepinfc*(1-lrcss)*r(-1) + (rbcss*grconstepinfc-rss)*lrcss/grconstepinfc*lrc(-1) ) - pinfc;
        z2c = nc - nc(-1) + z + alfac/(1-alfai)*v;

        // Definition of z2 for I sector
		//z2i = 1/((rbiss-rss)*lriss+rss) * ( rbiss*lriss*rbi + rss*(1-lriss)*r(-1) + (rbiss-rss)*lriss*lri(-1) );
        //z2i = grconstepinfc/((rbiss-rss)*lriss+rss) * ( rbiss*lriss*(rbi+pinfc) + rss/grconstepinfc*(1-lriss)*r(-1) + (rbiss*grconstepinfc-rss)*lriss/grconstepinfc*lri(-1) ) - pinfc;
        z2i = ni - ni(-1) + z + alfac/(1-alfai)*v;

        // Binding incentive constraint for bank in C sector
    	//qc + sc = lrc + nc;

        // Binding incentive constraint for bank in I sector
    	//qi + si = lri + ni;

        // Leverage ratio of bank in C sector
		lrc = etac + nucss/(clamb-nucss)*nuc;

        // Leverage ratio of bank in I sector
		lri = etai + nuiss/(clamb-nuiss)*nui;

        // Wealth accumulation equation of bank in C sector
    	//nc = sncss*cthetab*lrcss*exp(-ga-alfac/(1-alfai)*gv) * (rbcss*rbc + (1/lrcss-1)*rss*r(-1) 
        //    + (rbcss-rss)*lrc(-1)) +sncss*cthetab*exp(-ga-alfac/(1-alfai)*gv) * (  (rbcss-rss)*lrcss + rss  ) * (  -z-alfac/(1-alfai)*v+nc(-1)  ) 
        //    +(1-sncss*cthetab*exp(-ga-alfac/(1-alfai)*gv) * (  (rbcss-rss)*lrcss + rss ) ) * (qc+sc - psicss/(1-psicss)*psic ) 
        //    + snc;
        //(  cthetab*exp(-ga-alfac/(1-alfai)*gv)*((rbcss-rss)*lrcss+rss)  +( 1-cthetab*((rbcss-rss)*lrcss+rss)  ) ) * sn;

       // NEW - with stochastic cthetab:
       //nc = sncss*cthetabss*exp(-ga-alfac/(1-alfai)*gv)*((rbcss-rss)*lrcss+rss) * (nc(-1) + cthetab(-1) - z -alfac/(1-alfai)*v + snc)
       //     + sncss*cthetabss*exp(-ga-alfac/(1-alfai)*gv) * ( rbcss*lrcss*rbc - rss*lrcss*r(-1) + (rbcss-rss)*lrcss*lrc(-1) + rss*r(-1) )
       //     + (1-exp(-ga-alfac/(1-alfai)*gv)*((rbcss-rss)*lrcss+rss)*cthetabss*sncss) * (qc + sc + snc);
      // nc = sncss*cthetabss/grconstepinfc*exp(-ga-alfac/(1-alfai)*gv)*((rbcss*grconstepinfc-rss)*lrcss+rss) * (nc(-1) + cthetab(-1) - z -alfac/(1-alfai)*v - pinfc + snc)
      //      + sncss*cthetabss/grconstepinfc*exp(-ga-alfac/(1-alfai)*gv) * ( rbcss*lrcss*grconstepinfc*(rbc + pinfc) - rss*lrcss*r(-1) + (rbcss*grconstepinfc-rss)*lrcss*lrc(-1) + rss*r(-1) )
      //      + (1-exp(-ga-alfac/(1-alfai)*gv)*((rbcss*grconstepinfc-rss)*lrcss+rss)*cthetabss/grconstepinfc*sncss) * (qc + sc + snc);

   nc = sncss*cthetabss/grconstepinfc*exp(-ga-alfac/(1-alfai)*gv)*((rbcss*grconstepinfc-rss)*lrcss+rss) * (nc(-1) + cthetab(-1) - z -alfac/(1-alfai)*v - pinfc )
            + sncss*cthetabss/grconstepinfc*exp(-ga-alfac/(1-alfai)*gv) * ( rbcss*lrcss*grconstepinfc*(rbc + pinfc) - rss*lrcss*r(-1) + (rbcss*grconstepinfc-rss)*lrcss*lrc(-1) + rss*r(-1) )
            + (1-exp(-ga-alfac/(1-alfai)*gv)*((rbcss*grconstepinfc-rss)*lrcss+rss)*cthetabss/grconstepinfc*sncss) * (qc + sc ) +snc;
  
  // NOTE: normalized snc
        // Wealth accumulation equation of bank in I sector
    	//ni = sniss*cthetab*lriss*exp(-ga-alfac/(1-alfai)*gv) * (rbiss*rbi + (1/lriss-1)*rss*r(-1) + (rbiss-rss)*lri(-1))       
        //    +  sniss*cthetab*exp(-ga-alfac/(1-alfai)*gv) * (  (rbiss-rss)*lriss + rss  ) * (  -z-alfac/(1-alfai)*v+ni(-1)  )         
        //    + (1-sniss*cthetab*exp(-ga-alfac/(1-alfai)*gv) * (  (rbiss-rss)*lriss + rss ) ) * (qi+si - psiiss/(1-psiiss)*psii )  
        //    +  sni;

        //NEW - with stochastic cthetab:
        //ni = sniss*cthetabss*exp(-ga-alfac/(1-alfai)*gv)*((rbiss-rss)*lriss+rss) * (ni(-1) + cthetab(-1) - z -alfac/(1-alfai)*v + sni)
        //    + sniss*cthetabss*exp(-ga-alfac/(1-alfai)*gv) * ( rbiss*lriss*rbi - rss*lriss*r(-1) + (rbiss-rss)*lriss*lri(-1) + rss*r(-1) )
        //    + (1-exp(-ga-alfac/(1-alfai)*gv)*((rbiss-rss)*lriss+rss)*cthetabss*sniss) * (qi + si + sni);
        //ni = sniss*cthetabss/grconstepinfc*exp(-ga-alfac/(1-alfai)*gv)*((rbiss*grconstepinfc-rss)*lriss+rss) * (ni(-1) + cthetab(-1) - z -alfac/(1-alfai)*v - pinfc + sni)
        //    + sniss*cthetabss/grconstepinfc*exp(-ga-alfac/(1-alfai)*gv) * ( rbiss*grconstepinfc*lriss*(rbi  +pinfc) - rss*lriss*r(-1) + (rbiss*grconstepinfc-rss)*lriss*lri(-1) + rss*r(-1) )
        //    + (1-exp(-ga-alfac/(1-alfai)*gv)*((rbiss*grconstepinfc-rss)*lriss+rss)*cthetabss/grconstepinfc*sniss) * (qi + si + sni);

ni = sniss*cthetabss/grconstepinfc*exp(-ga-alfac/(1-alfai)*gv)*((rbiss*grconstepinfc-rss)*lriss+rss) * (ni(-1) + cthetab(-1) - z -alfac/(1-alfai)*v - pinfc )
            + sniss*cthetabss/grconstepinfc*exp(-ga-alfac/(1-alfai)*gv) * ( rbiss*grconstepinfc*lriss*(rbi  +pinfc) - rss*lriss*r(-1) + (rbiss*grconstepinfc-rss)*lriss*lri(-1) + rss*r(-1) )
            + (1-exp(-ga-alfac/(1-alfai)*gv)*((rbiss*grconstepinfc-rss)*lriss+rss)*cthetabss/grconstepinfc*sniss) * (qi + si ) +sni;

  // NOTE: normalized sni


//(  cthetab*exp(-ga-alfac/(1-alfai)*gv)*((rbiss-rss)*lriss+rss) + ( 1-cthetab*((rbiss-rss)*lriss+rss)  ) ) * sn;

        // Borrow in advance constraint for C sector
		kpc = sc;

        // Borrow in advance constraint for I sector
		kpi = si;

        // Leverage equation in C sector
        qc = lrc + nc - sc + psicss/(1-psicss)*psic;

        // Leverage equation in I sector
        qi = lri + ni - si + psiiss/(1-psiiss)*psii;

        // Stochastic return on assets for bank in C sector
    	rbc = 1/(crkc+qcss*(1-cdeltac)) * ( crkc*(rkc+zcapc) + qcss*(1-cdeltac)*(qc) ) -qc(-1) +spkc + z - (1-alfac)/(1-alfai)*v;

        rbcq = 1/(crkc+qcss*(1-cdeltac)) * (  qcss*(1-cdeltac)*(qc) ) -qc(-1) +spkc + z - (1-alfac)/(1-alfai)*v;

        // Stochastic return on assets for bank in I sector
    	rbi = 1/(crki+qiss*(1-cdeltai)) * ( crki*(rki+zcapi) + qiss*(1-cdeltai)*(qi) ) -qi(-1) +spki + z - (1-alfac)/(1-alfai)*v;

    	rbiq = 1/(crki+qiss*(1-cdeltai)) * (  qiss*(1-cdeltai)*(qi) ) -qi(-1) +spki + z - (1-alfac)/(1-alfai)*v;


        // Price of capital in C sector
		qc = exp(2*1/(1-alfai)*gv)*csadjcostc * ( invec-invec(-1)+1/(1-alfai)*v) - cbeta*exp(2*1/(1-alfai)*gv)*csadjcostc * (invec(+1)-invec+1/(1-alfai)*v(+1)) - qs + pi + (1+crho) * ( (iiss^(-crho)+icss^(-crho))^(-1) * (icss^(-crho)*invec + iiss^(-crho)*invei) - invec );

        // Price of capital in I sector
		qi = exp(2*1/(1-alfai)*gv)*csadjcostc * ( invei-invei(-1)+1/(1-alfai)*v) - cbeta*exp(2*1/(1-alfai)*gv)*csadjcostc * (invei(+1)-invei+1/(1-alfai)*v(+1)) - qs + pi + (1+crho) * ( (iiss^(-crho)+icss^(-crho))^(-1) * (icss^(-crho)*invec + iiss^(-crho)*invei) - invei );

        // External finance premium in C sector
        //spreadc = rbc(+1) - r;
		spreadc = rbcss*grconstepinfc/(rbcss*grconstepinfc - rss)*(rbc(+1) + pinfc(+1)) - rss/(rbcss*grconstepinfc - rss)*r;

        // External finance premium in I sector
        //spreadi = rbi(+1) - r;
		spreadi = rbiss*grconstepinfc/(rbiss*grconstepinfc - rss)*(rbi(+1) + pinfc(+1)) - rss/(rbiss*grconstepinfc - rss)*r;

        // Central Bank's feedback rule for C sector
        psic = cvarthetac*(rbc(+1)-r);

        // Central Bank's feedback rule for I sector
        psii = cvarthetai*(rbi(+1)-r);

        

    // Shock processes:
    // TFP growth shock to cons sector

        z = crhoz*z(-1)  + ez  + ez1aux(-1) + ez4aux(-4) + ez8aux(-8);
        ez1aux = ez1;
        ez4aux = ez4;
        ez8aux = ez8;
    
    // TFP stationary shock to cons sector
      
        z_l = crhozl*z_l(-1)  + ezl;
         
   // Investment-specific sector growth shock

        v = crhov*v(-1)  + ev + ev1aux(-1) + ev4aux(-4) + ev8aux(-8);
        ev1aux = ev1;
        ev4aux = ev4;
        ev8aux = ev8;      

    // Investment-specific stationary shock 
      
         v_l = crhovl*v_l(-1)  + evl;

   
   //  Preference, government, installation
   //  monetary policy, price markup, wage markup
	           
        b = crhob*b(-1) + eb;
	      
        g = crhog*g(-1) + eg;
        
        qs = crhoqs*qs(-1) + eqs;
	     
        ms = crhoms*ms(-1) + em;
	    
        spinfc = crhopinfc*spinfc(-1) + epinfmac - cmapc*epinfmac(-1);
	          epinfmac=epinfc;

        spinfi = crhopinfi*spinfi(-1) + epinfmai - cmapi*epinfmai(-1);
	          epinfmai=epinfi;
	    
        sw = crhow*sw(-1) + ewma - cmaw*ewma(-1) ;
	          ewma=ew; 
    // Shock to inflation objective
	      spi = crhospi*spi(-1) + espi;

    
    // Shock to bank's equity capital in C and I sector
       	snc = crhosnc*snc(-1) + enc;		
        sni = crhosni*sni(-1) + eni;
        //sn = crhosn*sn(-1) + en;	

    // Shock to the quality of physical capital in C and I sector
        spkc = crhospkc*spkc(-1) + epkc + epkc1aux(-1) + epkc2aux(-2) + epkc3aux(-3) + epkc4aux(-4) + epkc8aux(-8);
        epkc1aux = epkc1;
        epkc2aux = epkc2;
        epkc3aux = epkc3;
        epkc4aux = epkc4;
        epkc8aux = epkc8;
        spki = crhospki*spki(-1) + epki + epki1aux(-1) + epki2aux(-2) + epki3aux(-3) + epki4aux(-4) + epki8aux(-8);
        epki1aux = epki1;
        epki2aux = epki2;
        epki3aux = epki3;
        epki4aux = epki4;
        epki8aux = epki8;

    // Shock to bank's equity
       cthetab = crhothetab*cthetab(-1) + ethetab;

// TFP shock common to both sectors
        tfp = crhotfp*tfp(-1) + etfp + etfp4aux(-4) + etfp8aux(-8);
        etfp4aux = etfp4;
        etfp8aux = etfp8;


// measurement equations
/*
dy=y-y(-1) + z + alfac/(1-alfai)*v + (cga + alfac/(1-alfai)*cgv) ;
dc=c-c(-1) + z + alfac/(1-alfai)*v + (cga + alfac/(1-alfai)*cgv)  ;
dinve=inve-inve(-1) + 1/(1-alfai)*v + ( 1/(1-alfai)*cgv) ;
dpi = pi - pi(-1) + z + (alfac-1)/(1-alfai)*v + (cga+ (alfac-1)/(1-alfai)*cgv) ;
dw=w-w(-1) + z + alfac/(1-alfai)*v + (cga + alfac/(1-alfai)*cgv)  ;
pinfobsc = 1*(pinfc) + constepinfc ;
pinfobsi = 1*(pinfi) + constepinfi ;
robs =    r + log(constepinfc/100 + 1) - log(cbeta) + (cga + alfac/(1-alfai)*cgv) ;
labobs = lab + log(Lss) ;
rrate = r - pinfc(+1); 
spreadobsi = rbi(+1)+log(rbiss)-robs;
spreadobsc = rbc(+1)+log(rbcss)-robs;
dequity = exp(ga + alfac/(1-alfai)*gv) * ( ( (qcss*scss/lrcss)/(qcss*scss/lrcss+qiss*siss/lriss) ) *(nc-nc(-1)) 
+ (1- ( (qcss*scss/lrcss)/(qcss*scss/lrcss+qiss*siss/lriss) ) )*(ni-ni(-1)) + z + alfac/(1-alfai)*v +cga + alfac/(1-alfai)*cgv);
*/
// NEW measurement equations (using demeaned data)

dy=y-y(-1) + z + alfac/(1-alfai)*v   ;
dc=c-c(-1) + z + alfac/(1-alfai)*v   ;
dinve=inve-inve(-1) + 1/(1-alfai)*v  ;
dpi = pi - pi(-1) + z + (alfac-1)/(1-alfai)*v  ;
dw=w-w(-1) + z + alfac/(1-alfai)*v   ;
pinfobsc = 1*(pinfc)  ;
pinfobsi = 1*(pinfi)  ;
robs =    r  ;
labobs = lab + log(Lss) ;
rrate = r - pinfc(+1); 
spreadobsi = exp(rbi(+1))*exp(pinfc(+1)) - exp(robs);
spreadobsc = exp(rbc(+1))*exp(pinfc(+1)) - exp(robs);
dequity = exp(ga + alfac/(1-alfai)*gv) * ( ( (qcss*scss/lrcss)/(qcss*scss/lrcss+qiss*siss/lriss) ) *(nc-nc(-1)) 
+ (1- ( (qcss*scss/lrcss)/(qcss*scss/lrcss+qiss*siss/lriss) ) )*(ni-ni(-1)) + z + alfac/(1-alfai)*v );


// spread in the data = average spread in the model
// using qx*sx /nx = lrxdata ==> nx = qx*sx *lrxdata^-1
//spread = ( rbc(+1)+log(rbcss)-robs )* (( exp(qc+log(qcss))*exp(sc+log(scss)) )/exp(nc+log(qcss*scss*lrcdata^-1)) ) 
//        + ( rbi(+1)+log(rbiss)-robs )* (( exp(qi+log(qiss))*exp(si+log(siss)) )/exp(ni+log(qiss*siss*lridata^-1)) ) + esp1;
//dcredit = ( qc  - qc(-1) + sc - sc(-1) + z + alfac/(1-alfai)*v + cga + alfac/(1-alfai)*cgv )* ( ( exp(qc(-1)+log(qcss))*exp(sc(-1)+log(scss)) )/ (  exp(qc(-1)+log(qcss))*exp(sc(-1)+log(scss)) + exp(qi(-1)+log(qiss))*exp(si(-1)+log(siss)) ) )
//         + ( qi  - qi(-1) + si - si(-1) + z + alfac/(1-alfai)*v + cga + alfac/(1-alfai)*cgv )* ( ( exp(qi(-1)+log(qiss))*exp(si(-1)+log(siss)) )/ (  exp(qc(-1)+log(qcss))*exp(sc(-1)+log(scss)) + exp(qi(-1)+log(qiss))*exp(si(-1)+log(siss)) ) ) + esp2;

end; 

//steady;

shocks;

var ez;
stderr 0.5;
var ev;
stderr 0.5;
var eb;
stderr 0.1;
var eg;
stderr 0.5;
var eqs;
stderr 0.0;
var em;
stderr 0.1;
var epinfc;
stderr 0.1;
var epinfi;
stderr 0.1;
var ew;
stderr 0.1;
var espi;
stderr 0.0;
var enc;
stderr 0.0;
var eni;
stderr 0.0;
var ezl;
stderr 0.0;
var evl;
stderr 0.0;
var epki;
stderr 0.5;
var epkc;
stderr 0.5;
var epkc1;
stderr 0.0; 
var epkc2;
stderr 0.0; 
var epkc3;
stderr 0.0;
var epkc4; 
stderr 0.0;
var epkc8; 
stderr 0.0;
var epki1; 
stderr 0.0;
var epki2; 
stderr 0.0;
var epki3; 
stderr 0.0;
var epki4;
stderr 0.0;
var epki8;
stderr 0.0;
var ez1;
stderr 0.0;
var ez4;
stderr 0.0/sqrt(2);
var ez8;
stderr 0.0/sqrt(2);
var ev1;
stderr 0.0;
var ev4;
stderr 0.0/sqrt(2);
var ev8;
stderr 0.0/sqrt(2);
var ethetab;
stderr 0.0;
var etfp;
stderr 0.0;
var etfp4;
stderr 0.0/sqrt(2);
var etfp8;
stderr 0.0/sqrt(2);
end;




estimated_params;

// PARAM NAME, INITVAL, LB, UB, PRIOR_SHAPE, PRIOR_P1, PRIOR_P2, PRIOR_P3, PRIOR_P4, JSCALE
// PRIOR_SHAPE: BETA_PDF, GAMMA_PDF, NORMAL_PDF, INV_GAMMA_PDF




stderr ez,0.15,0.01,100,INV_GAMMA_PDF,0.5,2.0;
stderr ev,1.34,0.001,100,INV_GAMMA_PDF,0.5,2.0;
stderr eb,1.29,0.025,150,INV_GAMMA_PDF,0.1,2.0;
stderr eg,0.44,0.01,100,INV_GAMMA_PDF,0.5,2.0;
//stderr eqs,4.27,0.01,120,INV_GAMMA_PDF,0.5,2.0;
stderr em,0.11,0.01,80,INV_GAMMA_PDF,0.1,2.0;
//stderr espi,0.20,0.01,80,INV_GAMMA_PDF,0.1,2.0;
stderr epinfc,0.66,0.01,80,INV_GAMMA_PDF,0.1,2.0;
stderr epinfi,0.24,0.01,80,INV_GAMMA_PDF,0.1,2.0;
stderr ew,0.38,0.01,20,INV_GAMMA_PDF,0.1,2.0;
//stderr enc,0.260,0.01,100,INV_GAMMA_PDF,0.1,2.0;
//stderr eni,0.19,0.01,100,INV_GAMMA_PDF,0.1,2.0;
stderr epkc,0.1760,0.01,100,INV_GAMMA_PDF,0.5,2.0;
stderr epki,0.70,0.01,100,INV_GAMMA_PDF,0.5,2.0;
//stderr ethetab,0.001,0.0001,10,INV_GAMMA_PDF,0.001,2.0;
//stderr ezl,0.57,0.01,40,INV_GAMMA_PDF,0.5,2.0;
//stderr evl,0.95,0.01,40,INV_GAMMA_PDF,0.5,2.0;
//stderr esp1, 0.05,INV_GAMMA_PDF,0.05,0.01;
//stderr esp2, 0.05,INV_GAMMA_PDF,0.05,0.01;
//stderr epkc1,0.05,0.01,100,INV_GAMMA_PDF,0.5,2.0;
//stderr epkc2,0.05,0.01,100,INV_GAMMA_PDF,0.5,2.0;
//stderr epkc3,0.05,0.01,100,INV_GAMMA_PDF,0.5,2.0;
stderr epkc4,0.06,0.01,100,INV_GAMMA_PDF,0.5/sqrt(2),2.0;
stderr epkc8,0.205,0.01,100,INV_GAMMA_PDF,0.5/sqrt(2),2.0;
//stderr epki1,0.05,0.01,100,INV_GAMMA_PDF,0.5,2.0;
//stderr epki2,0.05,0.01,100,INV_GAMMA_PDF,0.5,2.0;
//stderr epki3,0.05,0.01,100,INV_GAMMA_PDF,0.5,2.0;
stderr epki4,0.12,0.01,100,INV_GAMMA_PDF,0.5/sqrt(2),2.0;
stderr epki8,0.15,0.01,100,INV_GAMMA_PDF,0.5/sqrt(2),2.0;



//stderr etfp,0.37,0.01,100,INV_GAMMA_PDF,0.5,2.0;

//stderr ez1,0.05,0.01,100,INV_GAMMA_PDF,0.1/sqrt(1),2.0;
//stderr ez4,0.415,0.01,100,INV_GAMMA_PDF,0.5/sqrt(2),2.0;
//stderr ez8,0.185,0.01,100,INV_GAMMA_PDF,0.5/sqrt(2),2.0;
//stderr ev1,0.05,0.01,100,INV_GAMMA_PDF,0.1/sqrt(1),2.0;
//stderr ev4,0.085,0.01,100,INV_GAMMA_PDF,0.5/sqrt(2),2.0;
//stderr ev8,0.49,0.01,100,INV_GAMMA_PDF,0.5/sqrt(2),2.0;

//stderr etfp4,0.33,0.01,100,INV_GAMMA_PDF,0.5/sqrt(2),2.0;
//stderr etfp8,0.12,0.01,100,INV_GAMMA_PDF,0.5/sqrt(2),2.0;


crhoz,0.78 ,.01,.9999,BETA_PDF,0.4,0.20;
crhov,.13 ,.01,.9999,BETA_PDF,0.4,0.20;
crhob,.91,.01,.9999,BETA_PDF,0.6,0.20;
crhog,.98,.01,.9999,BETA_PDF,0.6,0.20;
//crhoqs,.97,.01,.9999,BETA_PDF,0.6,0.20;
//crhoms,.38,.0000000001,.9999,BETA_PDF,0.4,0.20;
//crhospi,.95,.0000000001,.9999,BETA_PDF,0.4,0.20;
crhopinfc,.123,.01,.9999,BETA_PDF,0.6,0.20;
crhopinfi,.84,.01,.9999,BETA_PDF,0.6,0.20;
crhow,.0429,.001,.9999,BETA_PDF,0.6,0.20;
//crhosnc,0.82 ,.01,.9999,BETA_PDF,0.5,0.20;
//crhosni,0.55 ,.01,.9999,BETA_PDF,0.5,0.20;
crhospkc,0.82 ,.01,.9999,BETA_PDF,0.6,0.20;
crhospki,0.93 ,.01,.9999,BETA_PDF,0.6,0.20;
//crhozl,.936 ,.01,.9999,BETA_PDF,0.6,0.20;
//crhovl,.532 ,.01,.9999,BETA_PDF,0.6,0.20;
//crhothetab,0.7 ,.01,.9999,BETA_PDF,0.5,0.20;
//crhotfp,0.7 ,.01,.9999,BETA_PDF,0.6,0.20;


//cmapc,.28,0.01,.9999,BETA_PDF,0.5,0.2;
//cmapi,.60,0.01,.9999,BETA_PDF,0.5,0.2;
//cmaw,.95,0.01,.9999,BETA_PDF,0.5,0.2;


csadjcostc,1.46,0.2,50,GAMMA_PDF,4,1.0;
//csadjcosti,3.46,0.2,50,GAMMA_PDF,4,1.0;

chabb,0.567,0.001,0.99,BETA_PDF,0.5,0.1;
cprobw,0.71,0.01,0.99,BETA_PDF,0.66,0.1;
csigl,0.7,0.025,10,GAMMA_PDF,2,0.75;
cprobpc,0.70,0.01,0.99,BETA_PDF,0.66,0.10;
cprobpi,0.71,0.01,0.99,BETA_PDF,0.66,0.10;
cindw,0.10,0.000000001,0.99,BETA_PDF,0.5,0.15;
cindpc,0.1,0.000000001,0.99,BETA_PDF,0.5,0.15;
cindpi,0.23,0.000000001,0.99,BETA_PDF,0.5,0.15;
czcapi,4.892,0.000001,50,GAMMA_PDF,5.0,1.0;
czcapc,4.250,0.000001,50,GAMMA_PDF,5.0,1.0;
//czcap,12.2,0.000001,50,GAMMA_PDF,5.0,1.0;
//czcap,0.8,0.000001,1,BETA_PDF,0.5,0.1;

crpi,1.54,1.0,5,NORMAL_PDF,1.7,0.30;
crr,0.84,0.1,0.99,BETA_PDF,0.60,0.20;
//cry,0.09,0.0000000001,0.8,NORMAL_PDF,0.125,0.05;
crdy,0.75,0.000001,0.9,NORMAL_PDF,0.25,0.1;
//cry,0.16,0.0000000001,0.8,NORMAL_PDF,0.25,0.1;

//cminrho,0.1,0.00000000001,0.999999,BETA_PDF,0.5,0.20;



end;






//varobs dc dinve dy labobs dw  pinfobsc pinfobsi robs ; //spread dcredit;
varobs dc dinve dy labobs dw  pinfobsc pinfobsi robs spreadobsc spreadobsi dequity; //dequity spread dcredit;
// For Estimation
estimation(datafile= USdata_10, mode_compute=0,  mode_file= twosb_v12b_B_CQ48restrrho_mode, prior_trunc=1e-50, first_obs=1, presample=1, lik_init=1, prefilter=0, mh_replic=0,  mh_nblocks=2,mh_jscale=0.28,mh_drop=0.2 ) dc;
//estimation(datafile= USdata_6, mode_compute=4, mode_file = twosb2v6_mode, mode_check, prior_trunc=1e-50, first_obs=1, presample=1, lik_init=1, prefilter=0, mh_replic=0, mh_nblocks=2,mh_jscale=0.20,mh_drop=0.2, nodiagnostic) dc;
//estimation(datafile= USdata_6, mode_compute=0, mode_file = twosb2v6_mode, prior_trunc=1e-50, first_obs=1, presample=1, lik_init=1, prefilter=0, mh_replic=50000, mh_nblocks=2,mh_jscale=0.26,mh_drop=0.2) dc;
//load_mh_file,

// MDD after 50000:
// MDD 

//LogL = -746
//-719 with alpha = 0.36 and 1.15
//-724 with alpha = 0.3 and 1.15

// For IRFs
//estimation(datafile= USdata_6, mode_compute=0, mode_file = twosb2v6_mode, prior_trunc=1e-50, first_obs=1, presample=1, lik_init=1, prefilter=0, mh_replic=0, mh_nblocks=2,mh_jscale=0.25,mh_drop=0.2) dc;
//stoch_simul(irf = 40) dy dc dinve dw pinfobsc pinfobsi robs labobs spreadobsi spreadobsc;

//mode_file = twosb_v12b_B_CQ48restrrho_mode,
/*
cprobw=   0.01; 
clandawaux = 1.01-1; // ss wage mark up
cindw = 0.01;
cprobpc=   0.01; 
cprobpi = 0.01;
clandawaux = 1.01-1; // ss wage mark up
clandapaux = 1.01-1; // ss wage mark up

cindw = 0.01;
cindpc = 0.01;
cindpi = 0.01;

//crpi = 1.01;

clandap = 1 + clandapaux;
clandaw = 1 + clandawaux;
*/
//czcapc = 10.1;
//czcapi = 10.1;
//stoch_simul(conditional_variance_decomposition=[4 8 12 16 20 24 30 40], irf=0)  dc dinve dy labobs dw  pinfobsc pinfobsi robs spreadobsc spreadobsi dequity; 
//conditional_variance_decomposition=[1,4,8,12,16,20];
// conditional_variance_decomposition=[1 4 8]
//y inve c labobs labi labc spreadobsc spreadobsi robs pinfobsc dpi nc ni;